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We report results of systematic simulations of the dynamics of solitons in the framework of the one- 
dimensional nonlinear Schrodinger equation (NLSE), which includes the harmonic-oscillator (HO) potential 
and a random potential. The equation models experimentally relevant spatially disordered settings in Bose- 
Einstein condensates (BECs) and nonlinear optics. First, the generation of soliton arrays from a broad initial 
quasi-uniform state by the modulational instability (MI) is considered, following a sudden switch of the nonlin- 
earity from repulsive to attractive. Then, we study oscillations of a single soliton in this setting, which models a 
recently conducted experiment in BEC. Basic characteristics of the Mi-generated array, such as the number of 
solitons and their mobility, are reported as functions of the strength and correlation length of the disorder, and of 
the total norm. For the single oscillating soliton, its survival rate is found. Main features of these dependences 
are explained qualitatively. 
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I. INTRODUCTION 



The interplay of the disorder and nonlinearity is a topic which has been attracting a great deal of interest for a long time, see 
reviews Qb-Q]. In particular, much work has been done on the analysis of dynamics of solitons in disordered external potentials, 
chiefly in one-dimensional settings, see, e.g., Refs. H-d and references therein. These studies find an important application 
to the description of Bose-Einstein condensates (BECs) trapped in random potentials. The latter topic was addressed in many 
theoretical 0jl]-|[l3] and experimental |H]-||M] works. 

A ubiquitous theoretical model used in this field is the nonlinear Schrodinger equation, NLSE (alias the Gross-Pitaevskii 
equation, in terms of BEC If2llp for the wave function u(x, t), which includes a regular trapping potential and a term accounting 
for a disordered potential. The normalized form of this equation, with scaled time t and space coordinate x, is 

iu t + ^u xx + g\u\ 2 u - ]^VL 2 x 2 u - V dis (x)u = 0, (1) 

where g = +1 and —1 correspond to the self-attractive and self -repulsive nonlinearities, both signs being possible in BEC, 
characterizes the strength of the harmonic-oscillator (HO) trapping potential, and Vdi s (aO is a random potential representing 
the spatial disorder. In particular, this model describes dipole oscillations of the trapped condensate, experimentally studied in 
recent work 12011 . where strong effective dissipation induced by the disordered potential was discovered. On the other hand, 
Eq. (fl]i with t replaced by the propagation distance coordinate z is a model of a nonlinear optical waveguide with a random 
perturbation of the local refractive index, which is represented by the disordered potential |2||, hence the results reported below 
apply to spatial solitons in the nonlinear waveguide as well. 

Our objective is to systematically investigate fundamental properties of solitons in the framework of the model based on Eq. 
(H). These include the formation and motion of soliton trains due to the modulational instability (MI) of the broadly distributed 
condensate, and, as suggested by the experiment reported in Ref. 112011 . oscillations of a single soliton which is originally placed 
at a distance from the minimum of the trapping potential, x = 0. The results for these two types of the dynamical behavior, 
based on systematic simulations of Eq. (HJ and averaging over many realizations of the random setting, as well as on a qualitative 
analysis of the observed effects, are reported, respectively, in Sections II and III, and the paper is concluded by Section IV. 
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II. FORMATION OF SOLITON ARRAYS BY THE MODULATIONAL INSTABILITY IN THE SPATIALLY DISORDERED 

ENVIRONMENT 

The subject of this Section is the effect of the spatial disorder on the formation of soliton chains from an initial quasi-uniform 
state, and the subsequent evolution of the chains. The numerical simulations were subject to the following specifications. The 
spatial domain was taken as — 60ir < x < +607T, and the trapping strength is £1 2 = 25 x 10~ 4 , hence the respective HO length 
is much smaller than the size of the integration domain, which makes it possible to consider multi-soliton trains, and persistent 
motion of the solitons. The split-step Fourier method with a spatial grid composed of 1024 points and time step At = 0.001 was 
employed for the spatial and temporal discretization of Eg. 0}. The total integration time is T = 80. A fourth-order optimized 
implementation of the splitting was used, as in Refs. 122112311 . Further, the disorder was represented by a spatially correlated 
random function, V^ B (x) = Vdf(x,6), where Vd is the strength of the disorder, and the marginal distribution of f(x,9) has 
an exponential form, its covariance function being C(xi, x 2 ) — exp(—2(x 1 — x 2 ) 2 /V 2 ), in which V z is the correlation length. 
Other type of random functions as in M24I1 can be considered in a similar way. 

The initial quasi-uniform distribution of the condensate is taken as the ground-state solution of Eq. ([TJ with the self- repulsive 
nonlinearity, i.e., with g = —1, in the form of u — exp (—ifit) U(x), where /i is the real chemical potential, and function U(x) 
was found as a stationary solution of the associated nonlinear diffusion equation (the imaginary-time version J25ll of Eq. (Q])): 

Ut = -u xx - \u\ 2 u - ^Q 2 x 2 u - Vd is (x)u + fiu. (2) 

The same split-step method was used to solve Eq. (0. Then, the evolution of the MI (modulational instability) of this state 
was simulated in real time, following a sudden switch of the self -repulsion into self-attraction, i.e., replacement of g = — 1 by 
g = +1 in Eq. ((TJ, which exactly corresponds to the reversal of the sign of the nonlinearity by means of the Feshbach resonance 
in the well-known BEC experiment of 02611 . 

The MI splits the initial quasi-uniform state into a chain of solitons, which, generally speaking, is an obvious outcome of the 
evolution induced by the reversal of the nonlinearity sign. However, a new aspect of this outcome, on which we focus here, is 
the effect of the disordered environment on the resulting solitary wave chain. We have collected results which demonstrate the 
dependence of the characteristics of the emerging soliton chain on three control parameters: strength Vd and correlation length 
V z of the random potential, and the total norm of the initial state, N = J_°° \u(x)\ 2 dx. The characteristics whose variation was 
monitored are (i) the number of solitons; (ii) the largest displacement of solitons, i.e., the largest distance that peaks of individual 
solitons can travel in the course of the evolution; Basically the displacement for each soliton is computed as the distance between 
its leftmost and rightmost positions. Then the largest displacement is simply the largest value among all the solitons. (iii) the 
normalized average kinetic energy per soliton, which is defined as 

^=(E M i) _1 (^5 Mi ^ 2 )' (3) 

where the summation is performed on the full set of solitons in the emerging configuration, Vj is the velocity of the j-th soliton, 
and Mj — J \u(x)\ 2 dx is its effective mass, with the integration performed over a vicinity of the solitary wave's peak where the 
local amplitude of the field exceeds half of its peak value. 

The results presented below were produced by averaging over 100 different random realizations. Note that this leads to a 
non-integer number as the number of solitons for most cases. Typical examples of the realizations are displayed in Figs. Q]and 
|2j for small and large correlations lengths, V z = 5 and V z = 25, respectively. Individual solitons can be easily identified in the 
plots. 

The results for the number of solitons in emerging pattern and their largest displacement are summarized in Figs. [3] -[5] In 
each panel, we fix two of the above-mentioned control parameters (Vd, V z , N) and vary the third one. For example, the norm is 
varied in Fig. [3] each curve corresponding to a particular set of values of Vd and V z . 

The following conclusions can be made from Figs. [3]-[5] 

1. According to Fig. [3] the number of solitons in the chain increases linearly with the total norm, so that the mean norm 
per soliton, iV so i, is approximately constant. This is a direct effect of trapping the wave field by the random potential: in 
free space, the entire condensate tends to coalesce into a single soliton, which corresponds a minimum of the system's 
Hamiltonian ll27ll . However, the sufficiently strong disorder pins portions of the fragmented condensate, forcing them to 
self-trap into separated solitons. A fundamental threshold condition for the self-trapping of an initial state into an NLS 
soliton is known in the form of the condition imposed on the area of the initial configuration [27]: 

S= \u (x)\dx > So = In f 2 + V3J « 1.32. (4) 

J — oo 

For fixed parameters of the disorder, i.e., fixed average width L of local potential wells trapping fragments of the conden- 
sate, the above-mentioned constancy of the norm-per-soliton, N so \, implies a constant average amplitude, A ~ y N so \/L, 
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FIG. 1: (Color online) Two realizations of the random multi-soliton configuration generated by the development of the modulational instability 
for V z = 5. Here and in the next figure, top plots depict the spatiotemporal evolution of |w(a;,£)|, while the bottom ones display the 
corresponding total potential — the sum of the HO trap and the random part. 
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FIG. 2: (Color online) Two realizations of the random multi-soliton patterns for V z = 20 



hence the average area per soliton is constant too, S ~ AL ~ y N so \L, which is consistent with condition (|4j, that also 
implies an approximately constant area per soliton. Thus, condition (|4]i may explain the results observed in Fig. [3] 

2. Figure[3]demonstrates the increase of both the largest distance traveled by the solitons, and their mean kinetic energy, with 
the increase of the total norm. This feature may be explained by the fact that an individual soliton, moving through the 
disordered potential, is strongly braked due to the emission of radiation [1]. The effectively dissipative character of the 
motion of coherent wavepackets in this setting was recently demonstrated in the experiment of 1I20I1 (see also the following 
Section). However, if the system is filled by the trapped condensate, the moving soliton actually interacts with trapped 
segments of the condensate, i.e., with an effective pseudopotential [10], which is essentially smoother than the "bare" 
random potential. This effect leads to the suppression of the radiation losses, allowing the solitons to be more mobile. 

3. Figure|4]shows that the soliton number decreases, while the largest displacement increases, as the correlation length of the 
disorder, V z , increases. This conclusion is consistent with the conclusions presented in the previous item, as the increase 
of V z implies the transition to a less disordered potential. 

4. Figure|5]clearly shows that both the soliton number and largest displacement change very rapidly when Vd increases from 
zero to small finite values, i.e., the disorder starts to kick in when its strength is quite small. The observed jump of the 
soliton number to larger values, and the simultaneous drop of the free-path length are consistent with the above argument 



4 




15 20 25 30 35 40 45 15 20 25 30 35 40 45 

Norm Norm 



FIG. 3: (Color online) Results obtained for fixed strength and correlation length of the disorder, while the total norm varies. Here and in the 
two following figures, the left and right panels display, respectively, the largest displacement of solitons, and the number of solitons in the final 
configuration produced by the simulations. 
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FIG. 4: (Color online) The results for fixed total norm and strength of the disorder. 




FIG. 5: (Color online) The results for fixed total norm and correlation length of the disorder. 



5 



100 




15 20 25 30 35 40 45 
Norm 



FIG. 6: (Color online) The variation of the mean kinetic energy with increase of the total norm. 



which states that a deeper random potential splits the condensate into a large number of solitons, and impedes their free 
motion. 



III. OSCILLATIONS OF THE SOLITON IN THE COMBINED POTENTIAL 



Our next objective is to simulate the oscillatory motion of a single soliton in the combined HO an random potentials, which 
parallels the recently reported experiment performed in the condensate of 7 Li [20]. While the latter experiment focused on the 
repulsive interaction case of the (dissipative in the presence of disorder) dipolar motion of a full condensate, it is straightfor- 
ward to envision such dynamics for an attractive interaction condensate, namely a localized solitary wave. For this purpose, 
simulations of Eq. (Q]i were run with parameters selected as the rescaled version of those dealt with in the experiment, where 
the total number of atoms was ~ 10000, the scattering length a s is taken as three times the Bohr radius of Li, the transverse 
trapping frequency is uj r = 2tt x 260 Hz, and the longitudinal one is u z = 2ir x 5.5 Hz. The initial conditions were taken as 
uq(x) = -\/2aicsecri("\/2aic(;E — xq)), where xq is the initial shift of the soliton from the bottom of the HO potential, and a\c 
is determined by the total number of atoms. With the above physical values, f2 2 = 4.4749 x 10~ 4 and aic = 0.2269 were used 
in the simulations. The split-step Fourier method was employed in this case too, with a sufficiently small spacing in order to 
properly resolve the size of the soliton. In the course of the simulations, the disorder in Eq. (fTJ was turned on after one and a 
quarter of the period of the oscillations of the soliton in the HO potential, i.e., when the soliton's center was passing the origin, 
x = 0. 

For given initial shift xq, we mainly varied two parameters (as before), the correlation length V z and disorder strength Vd- For 
each fixed value of V z , we varied Vd to infer whether the soliton would survive after 10 periods of the oscillations. The goal 
was to compute the survival rate of the soliton under the influence of the disorder, using 10 realizations for this purpose. In 
each realization, the survival of the soliton was registered if its final norm, integrated over the FWHM range around its center, 
exceeded 50% of the initial value in the same range. The survival rate is then defined as the number of realizations featuring the 
surviving soliton, divided by 10 (the total number of the realizations). Typical examples of the survival and destruction of the 
oscillating solitons are displayed in Fig. [7] [Note that h in the horizontal axis label denotes the Planck constant]. 

The results for xq = 0.6 and 0.1 mm are displayed in Figs. [8] and [9] They correspond to the dimensionless values of 
xq = 254.6 and 42.4 in Eq.(Q]i. These clearly indicate that the survival rate drops to zero as the disorder strengthens and the 
correlation length decreases, in agreement with the experimental observations for the repulsive case ll20Tl . This is a natural 
consequence of the increasing rate of the emission of radiation by the soliton oscillating across the random potential. 

When the disorder becomes very strong, the soliton starts to survive again, as seen in Fig. [9] for xq = 0.1 mm. This is 
explained by the fact that a very strong random potential consists of local potential wells which quickly trap and immobilize the 
soliton, preventing its decay into radiation and, in addition, the strong random potential impedes the separation of the radiation 
waves from the parent soliton. The restabilization boundary is shown in Fig. [10] The approximately linear dependence of 
the minimum disorder strength, necessary for the restabilization, on Xo can be explained with the help the above argument: 
the driving force acting on the soliton in the HO potential grows linearly with .To, while the largest pinning force, induced by 
the random potential, is proportional to Vd- Therefore, the equilibrium between the two, which determines the restabilization 
threshold, implies Vd ~ xq . 
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FIG. 7: (Color online) The left and right panels display examples of solitons which, respectively, survive and get destroyed in the course 
of the oscillations, for an initial shift of the soliton to xo ~ 0.1 mm, its FWHM width 9.2142 /im, and a disorder correlation length of 
V z = 9.4264 /im. The disorder strength corresponding to the left and right panels is, respectively, V4 = h ■ 5(Hz) and h ■ 180(Hz), where h is 
the Planck constant. 
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FIG. 8: (Color online) The survival rate of the oscillating soliton for initial shift xq = 0.6 mm. 



IV. CONCLUSIONS 



The objective of this work was to report results of a systematic numerical analysis of individual soliton and soliton train dy- 
namics in the experimentally relevant setting based on the ID NLSE including the combination of the HO (harmonic-oscillator) 
and random potentials. This setting may be realized in BEC and nonlinear optics. Two basic problems were considered: the 
generation of soliton trains from the initial quasi-uniform state by the MI (modulational instability), after the sudden switch of 
the nonlinearity from the self-repulsion into attraction, and the survival or decay of a single soliton oscillating in this combined 
potential. Dependences of the main characteristics of the soliton train, and of the survival rate of the single oscillating soliton, on 
the strength and correlation length of the disordered potential, and also dependences of the characteristics of the soliton array on 
the total norm of the initial state, have been produced by averaging over a large number of random realizations. Salient features 
of this dependences have been explained in a qualitative form. 

A challenging problem is extending the analysis to the 2D setting, which may suggest new possibilities for the experiment. 
There, one has to consider the interplay of the above analyzed mechanisms with the potential collapse type events of super- 
critical atomic blobs and hence the relevant phenomenology will be considerable richer. The 2D setting is also of interest for 
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FIG. 9: (Color online) The survival rate of the oscillating soliton for xo = 0.1 mm. 



2500 



Threshold of the disorder strength 




0.1 0.2 0.3 0.4 0.5 0.6 



FIG. 10: (Color online) For given initial shift xo, the soliton survives if the disorder strength, Vj, exceeds the minimum value shown in this 
plot. 



the repulsive interaction dynamical case, whereby the formation of dark solitons and trains thereof reported in M20I1 may be 
substituted by the formation of vortices and vortex streets. 
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